/*
  File:    TAF.cpp
  Purpose: Utilities to evolve the tumor angiogenic growth factor. 

  Contact: Paul Macklin
           pmacklin@math.uci.edu
		  http://math.uci.edu/~pmacklin
*/

#include <cstdio>
#include <cstdlib>
#include <iostream>
#include <stdio.h>
#include <math.h>
using namespace std;

#include "../ImprovedMath.h"
#include "../Matrix.h"
#include "../DiffusionSolver.h"
#include "DiffusionSolverSteadyState.h"

Matrix TAF_Previous;

double OneFunctionBC( double )
{ return 1.0; }

double ZeroFunctionBC( double )
{ return 0.0; }

double TAFdiffusionCoefficientFunction( double taf,double x, double y,int i,int j)
{ 
 extern double TAFdiffusionConstant;
 extern Matrix Structure;
 
 if( *Structure(i,j) > 0.1 )
 { return 0.001; }  
 return TAFdiffusionConstant; 
}

double TAFsourceFunction(double taf, double x, double y,int i,int j)
{
 extern double _EPS;
 
 extern Matrix PhiQuiescent;
 extern Matrix PhiNecrotic;
 
 extern double TAFprodctionConstant;
 extern double TAFdecayConstant;

 double output = 0.0;

 if( *PhiQuiescent(i,j) <= _EPS && *PhiNecrotic(i,j) > _EPS  )
 { output += TAFprodctionConstant; }
 output -= TAFdecayConstant*taf;

 return output;
}

// we only degrade it and produce it in the hypoxic region. 
// The angio code produces it by the endothelial cells.

bool UpdateGlobalTAF_old( void )
{
 extern double dx;
 extern double dt;
 
 extern Matrix X;
 extern Matrix Y;
 extern int Xnodes;
 extern int Ynodes;
 
 extern Matrix TAF;

 extern double TAFdiffusionConstant;
 extern double TAFprodctionConstant;

 double s = 0;
 int factor = 1;
 double ds = dt/factor;
 while( ds > square(dx) )
 { factor++; ds = dt/factor; }
 
 double Tolerance = IntPow(dx,5); // IntPow(dx,4);
 double scale = MIN(TAFdiffusionConstant,TAFprodctionConstant);
 if( scale < 1 )
 { Tolerance *= scale; }

 for( int count = 0 ; count < factor ; count++ )
 {
  Matrix DiffusionConstant = ones(Xnodes,Ynodes);
  TAF = NonlinearDiffusionIteration( 
   X,Y,
   TAF,ds, 
   ZeroFunctionBC, 'E',  
   ZeroFunctionBC, 'E',     
   ZeroFunctionBC, 'E',   
   ZeroFunctionBC, 'E',  
   TAFdiffusionCoefficientFunction, 
   TAFsourceFunction,
   10000, Tolerance);
  s += ds;
 }
 return true;
}

// source functions for the quasi-steady solver

double TAFsourceTimesTAF( double taf, double x, double y, int i, int j )
{
 extern Matrix BloodVesselIndicator;
 extern Matrix SproutTips;
 extern Matrix VesselCells;
 extern Matrix CellColor;
 extern Matrix Oxygen;
 extern Matrix VEGFR;
 extern Matrix PhiQuiescent;
 extern Matrix PhiNecrotic;
 extern double TAFprodctionConstant;
 extern bool TumorBegin;
 extern double QuiescenceOxygenLevel;
 extern double TAFdecayConstant;
 extern double EndothelialTAFuptakeRate;
 
 extern double dx, dy;

 double output = 0.0;
 extern double _EPS; 
 // decay TAF everywhere
 
 output -= TAFdecayConstant;



if(TumorBegin)
{

 if( *PhiQuiescent(i,j) <= _EPS && *PhiNecrotic(i,j) > _EPS  )
 {
  output -= TAFprodctionConstant;
 }
 // note: all these results will be multiplied by taf inside the routine
  if( *SproutTips(i,j) > 0.25 )
 {
  output -= EndothelialTAFuptakeRate/(dx*dy);
 }
}
else
{
if( *BloodVesselIndicator(i,j) > 0.25&&(*CellColor(i,j)>15) )
 {
  output -= EndothelialTAFuptakeRate;
 }
 if(*Oxygen(i,j)<QuiescenceOxygenLevel)
 {
  output -= TAFprodctionConstant;
 }
}
 
 return output;
 
}

double TAFsourceRemainder( double taf, double x, double y, int i, int j )
{
 extern double _EPS;

 extern Matrix PhiQuiescent;
 extern Matrix PhiNecrotic;
  extern Matrix BloodVesselIndicator;
 extern Matrix SproutTips;
 extern Matrix CellColor;
 extern bool TumorBegin;
 extern Matrix Oxygen;
 extern double TAFprodctionConstant;
 extern double EndothelialTAFuptakeRate;
 double output = 0.0;
 extern double QuiescenceOxygenLevel;

 //if(i%5==0||j%5==0)
 //{
 //output = 1;
 //}

 // produce TAF in the hypoxic zone

// if( *PhiQuiescent(i,j) <= _EPS && *PhiNecrotic(i,j) > _EPS  )
// { output +=  } // TAFprodctionConstant; }

//if(*Oxygen(i,j)<0.5)
if(TumorBegin)
{
 if( *PhiQuiescent(i,j) <= _EPS && *PhiNecrotic(i,j) > _EPS  )
{
output += TAFprodctionConstant;
}
}
else
{
if(*Oxygen(i,j)<QuiescenceOxygenLevel)
{
output += TAFprodctionConstant;
}

}
//if(*BloodVesselIndicator(i,j)>0&&*CellColor(i,j)>15)
//{
//output += EndothelialTAFuptakeRate;
//}
 
 return output;
}

// the new, quasi-steady version

bool UpdateGlobalTAF( void )
{
 extern int Xnodes;
 extern int Ynodes;
 extern double dx;

 extern Matrix X;
 extern Matrix Y;
 
 extern Matrix TAF;
 
 extern Matrix PhiQuiescent;
 extern Matrix PhiNecrotic;
 extern Matrix Oxygen;
 extern double QuiescenceOxygenLevel;
 
 extern double _EPS;
 if( TAF_Previous.TellRows() == Xnodes && 
     TAF_Previous.TellCols() == Ynodes )
 {
  TAF = TAF_Previous;
 }
 
 bool IsHypoxic = false;
 for( int i=0 ; i < Xnodes ; i++ )
 {
  for( int j=0; j < Ynodes ; j++ )
  {
   if( *PhiQuiescent(i,j) <= _EPS && *PhiNecrotic(i,j) > _EPS )
   { IsHypoxic = true; }
   if( *Oxygen(i,j) < QuiescenceOxygenLevel && *PhiNecrotic(i,j) > _EPS )
   { IsHypoxic = true; }
  }
 }

 if( !IsHypoxic )
 {
  TAF = zeros(Xnodes,Ynodes);
  return true;
 }


 
 
 char BCtype = 'N';

 double v = (pow(5, 0.5) - 0.1)/(pow(5,0.5) - 1);
 
 /*
  for(int i=0; i< Xnodes; i++)
 {
 for(int j=0; j< Ynodes; j++)
 {
 //double r = pow( pow((j-50)*1.0/200,2) + pow((i-100)*1.0/200, 2) , 0.5);
 double r =  pow( pow( (i-(Xnodes-1)*1.0/2)/((Xnodes-1)*1.0/2),2) +pow( (j-(Ynodes-1)*1.0/2)/((Ynodes-1)*1.0/2),2) , 0.5);
 double a = -0.05;
 if(r>= 1)
 {
 *TAF(i, j) = 1;
 }
 else
 {
 *TAF(i, j) =exp(10*(-pow(r,3)/3+(a+1)*0.5*pow(r,2)-a*r)-10*(-1.0/3+(a+1)*0.5-a)); //exp(10*(-x.^3./3+(a+R).*0.5.*x.^2-a.*R.*x)-10*(-R.^3./3+(a+R).*0.5.*R.^2-a.*R.*R));

 }
  *TAF(i, j) = 0.001+1.0/200* j;//exp(-pow((j)*1.0/200,2)/0.45); 
 //cout<<*TAF(i, j) <<endl;

 }

 
 }
 */
 int MaxIterations;
 double Tolerance;


 MaxIterations = 10000; Tolerance = IntPow(dx,4); // 3);
 TAF = NonlinearDiffusionSteadyState(X,Y,TAF, 
    ZeroFunctionBC, 'N', ZeroFunctionBC, 'N', 
    ZeroFunctionBC, 'N', ZeroFunctionBC, 'N',  
    TAFdiffusionCoefficientFunction,
    TAFsourceTimesTAF,
    TAFsourceRemainder,
    MaxIterations, Tolerance );

  MaxIterations = 10000; Tolerance = IntPow(dx,5); // 4);
  TAF = NonlinearDiffusionSteadyState(X,Y,TAF, 
    ZeroFunctionBC, 'N', ZeroFunctionBC, 'N', 
    ZeroFunctionBC, 'N', ZeroFunctionBC, 'N',  
    TAFdiffusionCoefficientFunction,
    TAFsourceTimesTAF,
    TAFsourceRemainder,
    MaxIterations, Tolerance );
	    
 	
 TAF_Previous = TAF;
 


 return true;
}
